Trajectory optimization method and device for accurately deploying marine sensors under water

ABSTRACT

The disclosure discloses a trajectory optimization method and device for accurately deploying marine sensors under water. The method includes the following steps. 1. Randomly select N sets of initial control variables within a range. 2. Input all of N sets of x i  to the sensor&#39;s underwater glide kinematics and dynamics models, and calculate the smallest distance between N actual deployment positions and target deployment positions. 3. Determine whether the number of iterative operations is less than the preset value, if yes, perform global random walk and local random walk on N sets of x i , obtain N sets of x i  again, and return to step 2; otherwise, go to step 4. 4. Output the control variable x i  corresponding to Δs(x) nminmin  and the corresponding trajectory as the optimal control variable and optimal trajectory. The disclosure can improve the accuracy of prediction on the underwater three-dimensional trajectory of the marine sensor.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the priority benefit of China application Ser. No. 202110101658.1, filed on Jan. 26, 2021. The entirety of the above-mentioned patent application is hereby incorporated by reference herein and made a part of this specification.

BACKGROUND Technical Field

The disclosure belongs to the field of marine engineering and technology, and specifically relates to the field of deploying marine sensors under water, and more specifically, relates to a trajectory optimization method and device for accurately deploying marine sensors under water.

Description of Related Art

Most of the marine environment detection sensors are cylindrical or approximately cylindrical devices. They are commonly used in the field of marine monitoring and observation. They can measure and provide various marine environmental data. They are not only used in marine scientific research, but also serve as a crucial source of important data in the application of marine resource development. Taking the conductivity temperature depth (CTD) as an example, the CTD can be adapted to a variety of observation methods including fixed-point self-contained observation, disposable probe observation, etc., and obtain real-time high-precision physical environmental parameters such as temperature, raw data of conductivity/salinity and pressure/depth while carrying out investigation at the sea area water body scene. Therefore, CTD is currently one of the most commonly used instruments in physical oceanographic surveys. However, due to the lack of stable and reliable underwater point-to-point deploying technology directed for the CTD, it is difficult to ensure the accuracy of the deployment position of the CTD and the field observation range.

After the marine environment detection sensor is deployed, during the underwater gliding motion, the sensor may rotate around the axis itself and be affected by water. The motion and dynamic model for the marine environment detection sensor can be established to accurately predict the underwater motion trajectory of the marine environment detection sensor, which is also a basis for trajectory optimization for carrying out the deployment operation.

The paper “Motion Dynamics of Three-dimensional Bodies Falling Through Water” (Conference Paper “17th IWWWFB”, 2002, 21-81) discloses a three-dimensional motion model of falling objects without axial rotation. In the paper, external forces such as additional mass force and viscous force are taken into consideration in the three-dimensional space, but the lifting force that may be caused by the cylindrical axial rotation or the asymmetric vortex shedding is ignored. The paper “Numerical study on the trajectory of dropped cylindrical objects” (Ocean Engineering, 2017, 130: 1-9) discloses a three-dimensional motion model that takes into consideration of axial rotation of a cylinder, and the paper is used to simulate the underwater gliding motion trajectory of a cylinder. However, the lateral movement caused by asymmetric vortex shedding is ignored, and when dealing with the lifting force generated by the axial rotation of the cylinder, the ideal flow assumption is adopted to estimate based on the potential flow theory.

Neither of the two three-dimensional motion models calculates the lateral force caused by asymmetric vortex shedding. In particular, when the flow velocity of a cylinder is large, the lateral force caused by the asymmetric vortex is very large and the lateral motion will be very significant. As a result, the lateral dimension information is lost during trajectory prediction, resulting in huge prediction errors. Moreover, in the second file, the ideal flow assumption is adopted when calculating the lifting force generated by the axial rotation of the cylinder, ignoring the viscous effect of water, which makes the estimated lifting force larger than the actual value.

In addition, the above two comparative files only predict the trajectory of the cylinder without optimizing the trajectory. Therefore, the quality and accuracy of the deployment of marine environment detection sensor are not guaranteed when the task is performed under water.

SUMMARY

Based on the above shortcomings or needs of improvement for the related art, the disclosure completes high-precision point-to-point underwater deployment operations for marine sensors, and provides a trajectory optimization method and device for accurately deploying marine sensors under water, which conduct random update globally and random update locally directed at multiple input control variables. In the meantime, the multiple sets of multiple input control variables are simultaneously updated and compared in every iterative operation, thereby realizing trajectory optimization. The purpose of the disclosure is to improve the accuracy of three-dimensional trajectory prediction for underwater marine sensors, achieve trajectory optimization for specific point-to-point deployment operations, thereby realizing underwater marine sensor deployment task with high accuracy.

In order to achieve the above purpose, according to one aspect of the disclosure, a trajectory optimization method for accurately deploying a marine sensor under water is provided, which includes the following steps:

Step 1: According to the preset search range x_(min)≤x≤x_(max), N sets of values are randomly selected for the initial control variables x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(z0), ϕ₀, θ₀, ψ₀), denoted as x_(i). Specifically, i=1,2 . . . , N, x_(min), x_(max) are the lower limit and upper limit of the search range. U₀, V₀ and W₀ are the components of the initial linear velocity in the x-, y-, z-directions in the local coordinate system. ω_(x0), ω_(y0) and ω_(z0) are the components of the initial angular velocity ω₀ in the x-, y-, z-directions in the local coordinate system. (ϕ₀, θ₀, ψ₀) are the initial Euler angles of the sensor.

Step 2: In the n-th iterative operation, all N sets of x_(i) are input into the sensor's underwater glide kinematics and dynamics model, and N actual deployment positions are calculated. Then, the internal minimum deployment error Δs(x_(i))_(min) between target deployment position and the actual deployment position of each set of x_(i) in the current iterative operation is calculated. Thereafter, the minimum distances corresponding to the N sets of x_(i) are compared with each other. The minimum value is taken as the global minimum deployment error Δs(x)_(nminmin) of this iterative operation, wherein:

When n=1, Δs(x_(i))_(1min)=Δs (x_(i))₁, Δs (x)_(1minmin) is denoted as follows, and turn to step 3:

Δs(x)_(1minmin)=Min{Δs(x ₁)₁, Δs(x ₂)₁, . . . Δs(x _(N))₁};

When n>1, Δs (x_(i))_(nmin)=Min{Δs (x_(i))_(n), Δs(x_(i), Δs(x_(i))_((n−1)min)} is taken, and then turn to step 3. Under the circumstances,

Δs(x)_(nminmin)=Min{Δs (x ₁)_(nmin), Δs(x ₂)_(nmin), . . . Δs(x _(N))_(nmin)}

Specifically, the subscript n represents the current iteration operation steps, Δs(x_(i))_(n) represents the error distance between the actual deployment position and the target deployment position in the n-th step of iterative operation in the i-th set. The subscript min is the mark for minimum value in the set, and the subscript minmin is the mark for minimum value between the sets. Δs (x_(i))_(nmin) represents the internal minimum deployment error in the n-th step of iterative operation in the i-th set, and Δs(x)_(nminmin) represents the global minimum deployment error in the n-th step of iterative operation.

Step 3: It is determined whether n is less than the preset maximum number of iteration steps N₀, if yes, a global random walk is performed, N sets of x_(i) are updated, and then a local random walk is performed on the above updated N sets of x_(i) according to the preset probability. The N sets of x_(i) are obtained again, and return to step 2; otherwise, go to step 4.

Step 4: The control variable x corresponding to Δs(x)_(nminmin) and the corresponding trajectory are output as the optimal control variable and optimal trajectory.

Another technical problem to be solved by the disclosure mainly includes the inaccurate motion prediction on the lateral motion when the marine sensor performs underwater gliding motion along the axial direction, and the inaccurate motion prediction is firstly caused by the loss of lateral dimension motion information when the marine sensor performs underwater gliding motion without rotating in the axial direction. Secondly, when the marine sensor performs underwater gliding motion with rotating in the axial direction the lifting force value obtained from the inaccurate motion prediction is greater than the actual value, and the accuracy of deployment position and operation range is low.

To achieve the above purpose, furthermore, in step 2, the kinematics model of the sensor's underwater gliding motion is as follows.

$\overset{.}{U} = {{\frac{1}{m}\left\lbrack {{\left( {G - B} \right)\sin\;\theta} + F_{dx} + F_{ix}} \right\rbrack} - {W\;\omega_{y}} + {V\;\omega_{z}} + {L_{G}\left( {\omega_{y}^{2} + \omega_{z}^{2}} \right)}}$ $\overset{.}{V} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\;\theta\;\sin\;\phi} + F_{dy} + F_{ly} + F_{ey}} \right\rbrack} - {U\;\omega_{z}} - {L_{G}{\overset{.}{\omega}}_{z}}}$ $\overset{.}{W} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\;\theta\;\cos\;\phi} + F_{dz} + F_{lz} + F_{ez}} \right\rbrack} + {U\;\omega_{y}} + {L_{G}{\overset{.}{\omega}}_{y}}}$ ${\overset{.}{\omega}}_{y} = {\frac{1}{J_{y}}\left( {M_{ly} + M_{dy} + M_{ey} + M_{GBy} + {m\; L_{G}\overset{.}{W}} - {m\; L_{G}U\;\omega_{y}}} \right)}$ ${\overset{.}{\omega}}_{z} = {\frac{1}{J_{z}}\left( {M_{iz} + M_{dz} + M_{ez} - {m\; L_{G}\overset{.}{V}} - {m\; L_{G}U\;\omega_{z}}} \right)}$

In the formula, {dot over (U)}, {dot over (V)} and {dot over (W)} are the acceleration components of the sensor in the x-, y-, z-directions in the local coordinate system; {dot over (ω)}_(y) and {dot over (ω)}_(z) are the angular acceleration components of the sensor in the y-, z-directions in the local coordinate system.

m is the mass of the sensor, U, V, and W are the velocity components of the sensor in the x-, y-, z-directions in the local coordinate system. ω_(x), ω_(y) and ω_(z) are the angular acceleration components of the sensor in the x-, y-, z-directions in the local coordinate system. ϕ, θ, and ψ are the roll, pitch and yaw angles of the sensor in the local coordinate system. J_(y) and J_(z) are the moment of inertia of the sensor with respect to the y-axis and z-axis, respectively. B is the buoyancy force of the sensor. M_(GBy) is the hydrostatic moment caused by the relative offset L_(G) between the buoyancy center B and the gravity center G. F_(e) and M_(e) are the hydrodynamic force and moment of the sensor. F_(l) and M_(l) are the lifting force and torque received by the sensor. F_(d) and M_(d) are the drag force and torque received by the sensor. The subscripts x, y, z indicate that the corresponding parameter is the component on the corresponding coordinate axis.

There is a relationship between the angular velocity obtained through the conversion rule for the local coordinate system and the global coordinate system based on the Euler angle and the Euler angle:

$\begin{bmatrix} \omega_{x} \\ \omega_{y} \\ \omega_{z} \end{bmatrix} = {\begin{bmatrix} {{- \overset{.}{\psi}}\;{\sin(\theta)}} \\ {\overset{.}{\psi}\;{\sin(\phi)}\;\cos\;(\theta)} \\ {\overset{.}{\psi}\;{\cos(\phi)}\;{\cos(\theta)}} \end{bmatrix} + \begin{bmatrix} 0 \\ {\overset{.}{\theta}\;{\cos(\phi)}} \\ {{- \overset{.}{\theta}}\;\sin\;(\phi)} \end{bmatrix} + \begin{bmatrix} \overset{.}{\phi} \\ 0 \\ 0 \end{bmatrix}}$

Specifically, {dot over (ψ)} , {dot over (θ)} and {dot over (ϕ)} are the change rate of the yaw angle ψ, the pitch angle θ and the roll angle ϕ that varies with time.

Furthermore, the hydrodynamic forces in the y and z directions are calculated by the following expressions:

F _(ey) ={dot over (V)}m ₂₂−{dot over (ω)}_(y) m ₂₄ −UVM ₂₂(x _(t))+Uω _(z) x _(t) M ₂₂(x _(t))

F _(ez) ={dot over (W)}m ₃₃−{dot over (ω)}_(z) m ₃₅ −UVM ₃₃(x _(t))+Uω _(y) x _(t) M ₃₃(x _(t))

The hydrodynamic torques in the y and z directions are calculated according to the following formula:

M _(ey) =−{dot over (W)}m ₃₅−{dot over (ω)}_(y) m ₅₅ −UW[m ₃₃ +x _(t) mM ₃₃(x _(t))]+Uω _(y) [m ₃₅−(x _(t))² M ₃₃(x _(t))]

M _(ez) =−{dot over (V)}m ₂₆−{dot over (ω)}_(z) m ₆₆ −UV[m ₂₂ +x _(t) M ₂₂(x _(t))]+Uω _(z) [m ₂₆−(x _(t))² M ₂₂(x _(t))]

In the formula, m₂₂, m₂₄, m₂₆, m₃₃, m₃₅, m₅₅ and m₆₆ are the three-dimensional additional mass coefficients of the marine sensor, and M₂₂ (X_(t)) and M₃₃(x_(t)) respectively represent the two-dimensional additional mass coefficients of the effective trailing edge in the y and z directions. x_(t) is the longitudinal position of the effective trailing edge. For a sensor with a tip, M₃₃ (X_(t))=0, M₂₂ (x_(t))=0. For a marine sensor with two blunt ends, M₃₃(x_(t))≠0, M₂₂(x_(t))≠0.

Further, the drag force and torque are generated by the relative translational movement of the sensor in contact with water, and the axial drag force F_(dx) is calculated as follows:

$F_{dx} = \left\{ \begin{matrix} {{{{- C_{1}}\pi\; D\sqrt{v\;\rho^{2}L}U\sqrt{U}} - {\frac{1}{8}{\rho\pi}\; C_{dt}D^{2}U{U}}},} & {{\,_{\;}^{''}{Laminar}}\mspace{14mu}{flow}_{\;}^{''}} \\ {{{{- C_{2}}\frac{1}{2}{\rho\pi}\;{DLU}^{2}} - {\frac{1}{8}{\rho\pi}\; C_{dt}D^{2}U{U}}},} & {{\,_{\;}^{''}{Transition}}\mspace{14mu}{flow}_{\;}^{''}} \\ {{{C_{3}\frac{1}{2}{\rho\pi}\;{DLU}^{2}} - {\frac{1}{8}{\rho\pi}\; C_{dt}D^{2}U{U}}},} & {{}_{\;}^{}{}_{\;}^{}} \end{matrix} \right.$

In the expression, D is the cross-sectional diameter of the sensor, v is the kinematic viscosity of water, ρ is the density of water, C_(dt) is the tangential shape drag coefficient of the sensor, C₁, C₂, and C₃ are the shape drag coefficient of the sensor when there are laminar flow, transition flow and turbulence, L is the length of the sensor.

The lateral drag forces F_(dy) and F_(dz) in the y and z directions are calculated as follows:

F _(dy)=½ρC _(dy) D ∫ _(−0.5L) ^(0.5L)−(V+ω _(z) x)|V+ω _(z) x|dx

F _(dz)=½ρC _(dz) D ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)|W−ω _(y) x|dx

The resistance moments M_(dy) and M_(dz) caused by the lateral drag force are calculated as follows:

M _(dy)½ρC _(dz) D ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)|W−ω _(y) x|dx

M _(dz)½ρC _(dy) D ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)|V+ω _(z) x|dx

Specifically, C_(dy) and C_(dz) are respectively the drag coefficients in the y and z directions.

Further, if the marine sensor itself does not rotate around the axis, the lifting force F_(l) and the moment M_(l) are the lateral forces that are mainly caused by the asymmetry of the wake vortex distribution during the flow of water laterally crossing the sensor. If the sensor itself rotates, the water flow circulates laterally and generates pressure difference between the two sides of the sensor, which forms the lifting force.

The lifting force that is generated when the sensor does not rotate in the axial direction is caused by the asymmetry of the wake vortex distribution. The lifting forces are calculated as follows:

F _(ly)=½ρD ∫ _(−0.5L) ^(0.5L)(W−ω _(y) x)C _(l)(x)dx

F _(lz)=−½ρD ∫ _(−0.5L) ^(0.5L)(V+ω _(z) x)C _(l)(x)dx

M _(ly)=½ρD ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)C _(l)(x)dx

M _(lz)=½ρD ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)C _(l)(x)dx

Specifically, the lifting force coefficient C_(l)(x) is calculated as follows:

$\begin{matrix} {{C_{l}(x)} = \left\{ \begin{matrix} {{{KC}_{La}{\sin\left\lbrack {\pi\;{{St}\left( {t^{\prime} - 8} \right)}} \right\rbrack}},} & {t^{\prime} > 8} \\ {0,} & {t^{\prime} \leq 8} \end{matrix} \right.} & \; \end{matrix}$

In the expression, St=f_(v)D/V_(a) is the Strouhal number, f_(v) is the vortex shedding frequency, V_(a) is the cross-flow velocity, K=±1 reflects the random behavior of the lifting direction during gliding, C_(La) is the lift constant, and the time coefficient t′ is calculated as follows.

$\begin{matrix} {t^{\prime} = {\frac{s(x)}{R} - {{0.3}51}}} & \; \end{matrix}$

In the expression, s(x) is the displacement in the y direction when water flows around the cross section of the sensor, and R is the radius of the cross section of the sensor.

The lifting force and moment during axial rotation are calculated as follows:

F _(ly)=½ρDW ∫ _(−0.5L) ^(0.5L) C′ _(l)(x)(W−ω _(y) x)dx

F _(lz)=−½ρDV ∫ _(−0.5L) ^(0.5) C′ _(l) (x)(V+ω _(z) x)dx

M _(ly)=−½ρDV ∫ _(−0.5L) ^(0.5L) C′ _(l) 9 x)(V+ω _(x) x)xdx

M _(lz)=½ρDC′ _(l)(x)W ∫ _(−0.5L) ^(0.5L) (W−ω _(y) x)xdx

Specifically, the lifting force coefficient C′_(l)(x) is calculated as follows:

$\begin{matrix} {{{C_{l}^{\prime}(x)} = \left\{ \begin{matrix} {\frac{\omega_{x}D}{V_{a}},{\frac{\omega_{x}D}{V_{a}} \leq 8}} \\ {{8 + \ {{0.1}2\ \left( {\frac{\omega_{x}D}{V_{a}} - 8} \right)}}\ ,\ {\frac{\omega_{x}D}{V_{a}} > 8}} \end{matrix} \right.}.} & \; \end{matrix}$

Furthermore, the sensor in the local body coordinate system is an axisymmetric slender body, and the change of the initial roll angle will not affect the movement of the slender body. Therefore, the initial roll angle and the initial immersion depth are as follows:

ϕ₀=0

d=−L/2sin(θ₀)

Thus, the initial control variable is simplified as x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(z0), θ₀, ψ₀), and L represents the length of the sensor.

In order to achieve an efficient trajectory optimization algorithm for optimizing the target function, the algorithm has flexible iteration steps, which can be large or small, and thus being able to be free from the limitation of single step length, making the optimal solution search efficient. Further, in step 3, a global random walk is performed through Lévy flight, and global update is performed for N sets of control parameters x, as shown in the following formula:

x _(i) ^(n+1) =x _(i) ^(n) +a ⊗ Lévy(s, λ)

In the formula, a=a₀(x_(j) ^(t)−x_(i) ^(t)), a₀ is the preset constant, x_(j) represents the j-th set of control parameters that is randomly generated, s represents the step length of the iterative operation, ⊗ stands for item-based multiplication, and λ is the exponent of the Lévy distribution.

Further, after the global update in step 3, each set of control parameters x_(i) has the probability of p_(a)∈(0,1) to be locally updated. The update method of the local random walk is shown in the following formula:

x _(i) ^(n+1) =x _(i) ^(n+1)+ϵ·(x _(p) ^(n+1) −x _(q) ^(n+1)), p, q=1, 2, . . . N

Specifically, x_(p) and x_(q) are two random solutions, ϵ conforms to the uniform distribution of the interval [0,1], and the superscript n+1 is the number of iterations in the next step.

In order to achieve the above purpose, according to another aspect of the disclosure, a computer-readable/writable storage medium is provided, the computer-readable/writable storage medium stores a computer program. When the computer program is executed, the following steps are implemented including:

Step 1: randomly generate values for N sets of initial control variables such as velocity, accelerations and orientations within a certain search range. Input all of N sets of xi to the sensor's underwater glide kinematics and dynamics models, and calculate the smallest distance between N actual deployment positions and target deployment positions. 3. The disclosure can improve the accuracy of prediction on the underwater three-dimensional trajectory of the marine sensor.

Step 2: input the selected N sets of initial control variables into the sensor's underwater glide kinematics model.

Step 3: calculate the smallest distance between N actual deployment positions and target deployment positions.

Step 4: determine whether the number of iterative operations is larger than the pre-set value, if yes, output the trajectory and the corresponding control variable as the optimal trajectory and optimal control variable. Otherwise, perform global random walk and local random walk on N sets of control variable, obtain N sets of control variable again, and return to step 2.

Step 5: save the optimal trajectory and optimal control variable such as initial velocity, acceleration and orientation output in step 4 into the computer-readable/writable storage medium.

In order to achieve the above purpose, according to another aspect of the disclosure, a trajectory optimization device for accurately deploying a marine sensor under water is provided, which includes the computer-readable/writable storage medium and a processor and the computer-readable/writable storage medium is able to store the computer program and the data output by the execution of the computer program, the processor is used for calling and processing computer programs stored in the computer-readable storage medium. Once the computer program is called, it starts to solve the sensor's underwater glide kinematics model, output the underwater trajectory of the sensor and the deployment distance error. Finally, the optimal trajectory and optimal control variable such as initial velocity, acceleration and orientation can be obtained by sequential global and local random updates for control variable such as the initial velocity, acceleration and orientation.

In general, compared with the related art, the above technical solutions conceived in the disclosure can achieve the following advantageous effects:

1. The trajectory optimization algorithm designed by the disclosure greatly increases the number of main initial control variables, so that the target function can obtain a better optimal solution, and the distance between the deployment position and the actual target position is the smallest. In the meantime, the trajectory optimization algorithm of the disclosure expands the ability of trajectory optimization for multiple input variables. Through multiple sets of multiple input variables, the global random update and local random update are performed during each step of the iterative operation. In each step of the iterative operation, the prediction of the minimum error distance and comparison in a set and comparison between different sets are preformed, that is, multiple sets of multiple control variables are updated in parallel, comparison between different sets and within a set are performed to select the optimal result, so as to achieve a fast, efficient and accurate solution of the global optimal trajectory.

2. The underwater three-dimensional trajectory prediction model designed by the disclosure can take into account the asymmetric vortex shedding phenomenon that occurs along with the crossing flow, and calculate the lateral lifting force and lateral displacement of the marine sensor that does not rotate in axial rotation during the gliding process. Accordingly, the underwater three-dimensional trajectory prediction model provided by the disclosure increases the motion information of the lateral dimension, which is more in line with the actual three-dimensional gliding motion without axial rotation, and improves the engineering applicability of the prediction model.

3. The underwater three-dimensional trajectory prediction model designed by the disclosure modifies the calculation method of the lateral lifting force of the marine sensor during the gliding process while in axial rotation, and calculates lifting force generated by the marine sensor in viscous flow during axial rotation, and successfully takes into consideration the influence of the fluid viscosity. Therefore, the underwater three-dimensional trajectory prediction model provided by the present disclosure improves the accuracy of the motion information in the lateral dimension, and is more in line with the actual gliding motion of the marine sensor with axial rotation, which greatly improves the accuracy of the prediction model in engineering application.

4. The trajectory optimization algorithm designed by the disclosure has flexible iteration steps, which can be large or small, is able to be free from the limitation of single step length, and increases the search space per unit time. Through the global random walk algorithm and the local random walk algorithm, the main control parameters are updated twice, so that the main control parameters can be updated more frequently, which increases the efficiency of the optimal solution search.

5. The method for optimizing the underwater trajectory of the marine sensor designed by the disclosure can be applied to the trajectory optimization problems of the marine sensor in a wider range of operating requirements, including short time, high speed, short range, etc., and can be used to provide a more accurate guidance for optimal design and implementation of underwater deployment operation of marine sensors. In the meantime, the underwater trajectory optimization method designed by the disclosure can be used to solve the trajectory optimization problems for a variety of different types of marine sensors such as marine environment detection, underwater seismic monitoring, acoustic detection sensors, etc., and other underwater long-term observation type unpowered devices.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is the definition diagram of the coordinate system.

FIG. 2 is a motion information diagram of the marine environment detection sensor at lateral dimension when not in axial rotation during the gliding process.

FIG. 3 is a motion information diagram of the marine environment detection sensor at lateral dimension when in axial rotation during the gliding process.

FIG. 4 is a display diagram of the underwater deployment operation of the marine environment detection sensor.

FIG. 5 is a flow chart of the underwater trajectory optimization of the marine environment detection sensor.

FIG. 6 is the motion trajectory diagram of the marine environment detection sensor after optimization under the target function of the minimum deployment error.

DESCRIPTION OF THE EMBODIMENTS

In order to make the purpose, technical solutions, and advantages of the disclosure clearer, the following further describes the disclosure in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the disclosure, but not to limit the disclosure. In addition, the technical features involved in the various embodiments of the disclosure described below can be combined with each other as long as they do not conflict with each other.

As shown in FIG. 1 and FIG. 4, the marine environment detection sensor is deployed at a certain position underwater and performs a gliding motion until the sensor reaches a certain target position at the bottom of the water. In this embodiment, obtaining the minimum deployment error serves as an example, the underwater gliding motion trajectory of the marine environment detection sensor is optimized through the trajectory optimization method to set the main control variable value before the deployment. It should be noted that the trajectory optimization method is also applicable to the trajectory optimization of other marine sensors such as underwater seismic monitoring sensors, acoustic detection sensors, etc., and other underwater long-term observation type unpowered devices for underwater accurate deployment operations.

First, the underwater glide kinematics model is established first. As shown in FIG. 1, the global coordinate system used as the inertial system is {I} -OXYZ, and the coordinate system is regarded as positive when the Z axis is at the top. In the meantime, the local coordinate system is represented by {B}-oxyz, which is a right-hand Cartesian coordinate system fixed at the gravity center of the marine environment detection sensor.

In order to describe the motion of the sensor (in this embodiment, a slender cylinder) as a free fall in water, the following equations need to be solved in a local coordinate system.

$\overset{.}{U} = {{\frac{1}{m}\left\lbrack {{\left( {G - B} \right)\sin\theta} + F_{dx} + F_{ex}} \right\rbrack} - {W\omega_{y}} + {V\omega_{z}} + {L_{G}\left( {\omega_{y}^{2} + \omega_{z}^{2}} \right)}}$ $\overset{.}{V} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\theta\sin\phi} + F_{dy} + F_{ly} + F_{ey}} \right\rbrack} - {U\omega_{z}} - {L_{G}{\overset{.}{\omega}}_{z}}}$ $\overset{.}{W} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\theta\cos\phi} + F_{dz} + F_{lz} + F_{ez}} \right\rbrack} + {U\omega_{y}} + {L_{G}{\overset{.}{\omega}}_{y}}}$ ${\overset{.}{\omega}}_{y} = {\frac{1}{J_{y}}\left( {M_{ly} + M_{dy} + M_{ey} + M_{GBy} + {mL_{G}\overset{.}{W}} - {mL_{G}U\omega_{y}}} \right)}$ ${\overset{.}{\omega}}_{z} = {\frac{1}{J_{z}}\left( {M_{lz} + M_{dz} + M_{ez} - {mL_{G}\overset{.}{V}} - {mL_{G}U\omega_{z}}} \right)}$

In the formula, {dot over (U)}, {dot over (V)} and {dot over (W)} are the acceleration components of the sensor in the x-, y-, z-directions in the local coordinate system; {dot over (ω)}_(y) and {dot over (ω)}_(z) are the angular acceleration components of the sensor in the y-, z-directions in the local coordinate system. m is the mass of the marine environment detection sensor, U, V, and W are the velocity components of the sensor in translation in the x-, y-, z-directions in the local coordinate system. ω_(x), ω_(y) and ω_(z) are the angular velocity of the sensor in the x-, y-, z-directions in the local coordinate system. ϕ, θ, ψ are the Euler angles of the sensor in the local coordinate system, which represent the axial rotation angle, pitch angle, and yaw angle, respectively. J_(y) and J_(z) are the moment of inertia of the sensor with respect to the y-axis and z-axis, respectively. Buoyancy B is the hydrostatic pressure item related to the drainage volume of the marine environment detection sensor. M_(GBy) is the hydrostatic moment caused by the relative offset L_(G) between the buoyancy center B and the gravity center G. F_(e) and M_(e) are the hydrodynamic force and moment of the sensor. F_(l) and M_(l) are the lifting force and torque received by the sensor. F_(d) and M_(d) are the drag force and torque received by the sensor. The subscripts x, y, z indicate that the corresponding parameter is the component on the corresponding coordinate axis.

It should be noted that if the axial rotation rate along the x-axis is zero, that is, ω_(x)=0 at the initial moment, the axial rotation rate is always very close to zero during the falling process. Otherwise, suppose ω_(x)=c, and c is a fixed value throughout the process. In other words, in this embodiment, because the freely falling process is very short, for simplification, the possible deceleration of the axial rotation due to the viscosity of water is ignored.

The relationship between the angular velocity of the marine environment detection sensor and the Euler angle is obtained through the conversion rule of the local coordinate system and the global coordinate system as follows:

$\begin{bmatrix} \omega_{x} \\ \omega_{y} \\ \omega_{z} \end{bmatrix} = {\begin{bmatrix} {{- \overset{.}{\psi}}{\sin(\theta)}} \\ {\overset{.}{\psi}{\sin(\phi)}{\cos(\theta)}} \\ {\overset{.}{\psi}{\cos(\phi)}{\cos(\theta)}} \end{bmatrix} + \begin{bmatrix} 0 \\ {\overset{\overset{.}{\;}}{\theta}{\cos(\phi)}} \\ {{- \overset{\overset{.}{\;}}{\theta}}{\sin(\phi)}} \end{bmatrix} + \begin{bmatrix} \overset{.}{\phi} \\ 0 \\ 0 \end{bmatrix}}$

A dynamic model is established for each force in the underwater gliding motion equation, and the kinematics equation and dynamics model are solved at the same time to obtain each transient motion parameter and magnitude of force. The underwater glide kinematics model and dynamics model can be established offline and then directly called online during the trajectory optimization process.

(1) Calculation of hydrodynamic force and moment:

According to the potential flow theory under the slender body hypothesis, the hydrodynamic forces in the y and z directions are calculated as:

F _(ey) =−{dot over (V)}m ₂₂−{dot over (ω)} _(y) m ₂₄ −UVM ₂₂(x _(t))+Uω _(z) x _(t) M ₂₂(x _(t))

F _(ez) ={dot over (W)}m ₃₃−{dot over (ω)} _(z) m ₃₅ −UWM ₃₃(x _(t))+Uω _(y) x _(t) M ₃₃(x _(t))

The hydrodynamic torques in the y and z directions are calculated as follows:

M _(ey) ={dot over (W)}m ₃₅−{dot over (ω)} _(y)m₅₅ −UW[m ₃₃ +x _(t) mM ₃₃(x _(t))]+Uω _(y) [m ₃₅−(x _(t))² M ₃₃(x _(t))]

M _(ez) =−{dot over (V)}m ₂₆ −{dot over (ω)} _(z) m ₆₆ −UV[m ₂₂ +x _(t) M ₂₂ (x _(t))]+Uω _(z) [m ₂₆−(x _(t))² M ₂₂(x _(t))]

In the formula, m₂₂, m₂₄, m₂₆, m₃₃, m₃₅, m₅₅ and m₆₆ are the three-dimensional additional mass coefficients of the marine sensor, and M₂₂ (X_(t)) and M₃₃(x_(t)) respectively represent the two-dimensional additional mass coefficients of the effective trailing edge in the y and z directions. x_(t) is the longitudinal position of the effective trailing edge. For a marine environment detection sensor with a tip, the lateral force is zero: M₃₃ (x_(t))=0 M₂₂ (x_(t))=0. For a marine environment detection sensor with two blunt ends, M₃₃ (X_(t))≠0, M₂₂ (X_(t))≠0.

(2) Calculation of axial drag force: since the water flow is generated axially, and is characterized in that its contact area with the marine environment detection sensor is large for a long time, and is mainly composed of frictional resistance and shape resistance, the axial drag force F_(dx) is calculated as follows:

$F_{dx} = \left\{ \begin{matrix} {{{{- C_{1}}\pi D\sqrt{v\rho^{2}L}U\sqrt{|U|}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{{Laminar}\mspace{14mu}{flow}^{''}}} \\ {{{{- C_{2}}\frac{1}{2}\rho\pi DLU^{2}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{{Transition}\mspace{14mu}{flow}^{''}}} \\ {{{C_{3}\frac{1}{2}\rho\pi DLU^{2}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{Turbulence}^{''}} \end{matrix} \right.$

In the expression, D is the cross-sectional diameter of the marine environment detection sensor, v is the kinematic coefficient of water, ρ is the mass density of water, C_(dt) is the tangential shape drag coefficient of the sensor, C₁, C₂, and C₃ are the shape drag coefficient of the sensor when there are laminar flow, transition flow and turbulence, L is the length of the sensor.

(3) Lateral drag force: Since water flow is generated in the manner of laterally surrounding the marine environment detection sensor, and is characterized in that, since the sensor has a large aspect ratio, it is estimated by the strip theory. The lateral drag force F_(dy) and F_(dz) are calculated as follows:

F _(dy)=½ρC _(dy) D ∫ _(−0.5L) ^(0.5L)−(V+ω _(z) x)|V+ω _(z) x|dx

F _(dz)=½ρC _(dz) D ∫ _(−0.5L) ^(0.5L)−(W−ω _(y) x)|W−ω _(y) x)|W−ω _(y) x|dx

Correspondingly, the resistance moments M_(dy) and M_(dz) that is caused by lateral drag force have the following expressions:

M _(dy)=½ρC _(dx) D ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)|W−ω _(y) x|dx

M _(dz)=½ρC _(dy) D ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)|V+ω _(z) x|dx

In the formula, C_(dy) and C_(dz) are the drag coefficients in the y and z directions, respectively.

(4) If the marine environment detection sensor does not rotate axially, the lifting force and torque are the lateral forces that are generated mainly due to the wake vortex distribution asymmetry that is caused by the water flow that laterally surrounds the marine environment detection sensor. For the lifting force F_(ly) and the lifting force F_(lz) caused by the asymmetry of the wake vortex distribution, the lifting force moments M_(ly), and M_(lz) are calculated as follows:

F _(ly)=½ρD ∫ _(−0.5L) ^(0.5L) (W−ω _(y) x)C _(l)(x)dx

F _(lz)=−½ρD ∫ _(−0.5L) ^(0.5L) (V+ω _(z) x)C _(l)(x)dx

M _(ly)=½ρD ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)C _(l)(x)dx

M _(lz)=½ρD ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)C _(l)(x)dx

In the formula, the lifting force coefficient C₁(x) is calculated as follows:

$\begin{matrix} {{C_{l}(x)} = \left\{ \begin{matrix} {{KC_{La}{\sin\left\lbrack {\pi\;{{St}\left( {t^{\prime} - 8} \right)}} \right\rbrack}},} & {t^{\prime} > 8} \\ {0,} & {t^{\prime} \leq 8} \end{matrix} \right.} & \; \end{matrix}$

In the expression, St=f_(v)D/V_(a) is the Strouhal number, f_(v) is the vortex shedding frequency, V_(a) is the cross-flow velocity. It should be noted that K=±1 reflects the random behavior of the lifting direction during gliding. In the deployment that is performed at different initial pitch angles as shown in FIG. 2, the lateral movement trajectories of the marine environment detection sensor are shown in the form of vibrations. C_(La) is the lift constant, and the time coefficient t′ calculated is as follows.

$\begin{matrix} {t^{\prime} = {\frac{s(x)}{R} - {{0.3}51}}} & \; \end{matrix}$

In the formula, s(x) is the displacement in the y direction when water flows around the cross section of the sensor, and R is the radius of the cross section of the sensor.

(5) If the marine environment detection sensor rotates axially, the pressure difference between the two sides of the marine environment detection sensor is generated when the water flows laterally, which forms lifting forces F_(ly) and F_(lz), and lift moments are M_(ly) and M_(lz). FIG. 3 shows the changes in the three-dimensional motion trajectory of the marine environment detection sensor under different axial rotation rates.

F _(ly)=½ρDW ∫ _(−0.5L) ^(0.5L) C′ _(l) (x)(W−ω _(y) x)dx

F _(lz)−½ρDV ∫ _(−0.5L) ^(0.5L) C′ _(l) (x)(V+ω _(z) x)dx

M _(ly)=½ρDV ∫ _(−0.5L) ^(0.5L) C′ _(l)(x)(V+ω _(z) x)xdx

M _(lz)=½ρDC′ _(l)(x)W ∫ _(−0.5L) ^(0.5L) (W−ω _(y) x)xdx

In the formula, the lifting force coefficient C′_(l) (x) is calculated as follows:

$\begin{matrix} {{C_{l}^{\prime}(x)} = \left\{ \begin{matrix} {\frac{\omega_{x}D}{V_{a}},} & {\frac{\omega_{x}D}{V_{a}} \leq 8} \\ {{8 + {0.12\left( {\frac{\omega_{x}D}{v_{a}} - 8} \right)}},} & {\frac{\omega_{x}D}{V_{a}} > 8} \end{matrix} \right.} & \; \end{matrix}$

FIG. 5 shows the process of optimizing the underwater trajectory of the marine environment detection sensor. First, based on the target function value output by the trajectory prediction model under the N sets of main initial control variables, combined with the global and local random walking algorithms, the main initial control variables are updated, and iterative operation is performed in sequence until the number of updates reaches the preset maximum number of steps. The specific steps are as follows:

Step 1: The main initial control variables and initialization are determined. The initial control variables of the underwater three-dimensional trajectory of the marine environment detection sensor include 9 initial kinematics parameters x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(z0), ϕ₀, θ₀, ψ₀) and the initial position parameters (0, 0, d) of the marine environment detection sensor. U₀, V₀, W₀ are the components of the initial linear velocity in the three coordinate axis directions in the local coordinate system, ω_(x0), ω_(y0), ω_(z0) are the components of the initial angular velocity in the three coordinate axis directions in the local coordinate system, ϕ₀, θ₀, ψ₀ are the initial Euler angles, and d is the initial immersion depth.

In this embodiment, since the slender body is axisymmetric, the change of the initial roll angle will not affect the movement of the slender body, so the initial roll angle and the initial immersion depth are as follows:

ϕ₀=0

d=−L/2sin(θ₀)

To sum up, in the scenario specified in this embodiment, first it is determined that the final eight main initial control variables of the marine environment detection sensor include:

x=(U ₀ , V0, W ₀, ω_(x0), ω_(y0), ω_(z0), θ₀, ψ₀).

Then, for specific task requirements, the initial control variables should satisfy a certain search range, as shown in the following formula:

x _(min) ≤x≤x _(max)

N sets of x_(i) are selected randomly within the above range, wherein i=1,2 . . . , N, then the initialization is completed in the first step.

Step 2: All of N sets of x_(i) are input into the underwater glide kinematics model and dynamics model of the sensor, and N actual deployment positions are calculated. Then the minimum distance between the target deployment position and the actual deployment position within each set are calculated, and the minimum distances in the N sets are compared with each other, and the minimum value is denoted as:

Δs(x)_(nminmin)=Min{Δs(x ₁)_(nmin), Δs(x ₂)_(nmin), . . . Δs(x _(N))_(nmin})

In the formula, the subscript n represents the current iteration operation steps, Δs (x_(i))_(n) represents the error distance between the actual deployment position and the target deployment position in the n-th step of iterative operation in the i-th set. The subscript min is the mark for minimum value in the set, and the subscript minmin is the mark for minimum value between the sets. Δs (x_(i))_(nmin) represents the internal minimum deployment error in the n-th step of iterative operation in the i-th set, and Δs(x)_(nminmin) represents the global minimum deployment error in the n-th step of iterative operation.

When n=1, Δs(x_(i))₁ is directly used as Δs(x_(i))_(1min), Δs(x_(i))_(1minmin) is solved, and then go to step 3;

When n>1, Δs (x_(i))_(nmin)=Min{Δs (x_(i))_(n), Δs (x_(i))_((n−1)min)} is taken, then go to step 3.

Step 2: The target function of the smallest deployment errors are obtained and compared. The global minimum deployment error is the minimum distance between the actual deployment position and the target deployment position, as shown in the following formula:

Δs(x)_(nminmin) =Min{Δs(x ₁)_(nmin), Δs(x ₁)_(nmin) , Δs(x ₂)_(nmin), . . . Δs(x _(N))_(nmin)}

The N sets of x_(i) are input into the underwater gliding motion and dynamics model to get the deployment error Δs (x_(i))_(n) calculated in the n-th step. For step n=1, the initialized N sets of x_(i) are input into the underwater gliding motion and dynamics model to obtain the deployment error ds(x_(i))₁ calculated in step 1, and the obtained Δs(x_(i))₁ is retained and used as the initial deployment error Δs(x_(i))_(1min) in the i-th set, N sets of Δs(x_(i))_(1min) are compared, and the minimum value Δs(x)_(nminmin) is taken as the initial deployment error.

For step n>1, by comparing the Δs(x_(i))_(n) calculated in the n-step and the corresponding Δs(x_(i))_((n−1)min) in the n−1-th step, the minimum value is selected as the internal minimum deployment error Δs(x_(i))_(nmin) in the i-th set in the n-step. N sets of Δs(x_(i))_(nmin) are compared, and the minimum value Δs(x)_(nminmin) is taken as the global minimum deployment error in the n-th step, which is denoted as:

Δs(x)_(nminmin)=Min{Δs(x ₁)_(nmin) , Δs(x ₂)_(nmin) , . . . Δs(x _(N))_(nmin)}

Specifically,

Δs(x _(i))_(nmin)=Min{Δs(x _(i))_(n) , Δs(x_(i))_((n−1))min}

Δs(x _(i))_(n)=√{square root over ((X _(t)(x _(i))_(n) −a)²+(Y _(t)(x _(i))_(n) −b)²+(Z _(t)(x _(i))_(n) −c)²)}

(a,b,c) is the target deployment position, (X_(t)(x_(i))_(n), Y_(t)(x_(i))_(n), Z_(t)(x_(i))_(n)) is the actual deployment position at the n-th step of iterative operation in the i-th set.

Step 3: It is determined whether n is less than the preset maximum number of iteration steps N₀. If n<N₀, the N sets of initial control variables x in the n+1th step need to be updated through an optimization algorithm. The specific algorithm is as follows:

(3.1) Global random walk algorithm. A global random walk is performed through Lévy flight, and the first update is performed for the N sets of control parameters x_(i) , as shown in the following formula:

x _(i) ^(n+1) =x _(i) ^(n) +a ⊗ Lévy(s, λ),

In the formula, a=a₀(x_(j) ^(n)−x_(i) ^(n)), and a₀ is a constant. In most cases, a=O(1) can be adopted, that is, the order of a is usually 1. ⊗ stands for item-based multiplication, and λ is the exponent of the Lévy distribution. The advantage of using the Lévy flight to perform a global random walk in this embodiment is that its step length is flexible, which can be large or small, free from the limitation of a single step length, and increases the search space per unit time. In other embodiments, other common global update algorithms can also be used instead.

(3.2) Local random walk algorithm. After the first update (i.e., global update), there is a probability of p_(a) ∈ (0, 1) for each set of control parameters x_(i) in the N sets to be updated for the second time, mainly through a local random walk, as shown in the following formula:

x _(i) ^(n+1) =x _(i) ^(n+1)+ϵ·(x _(p) ^(n+1) −x _(q) ^(n+1)), p,q=1,2, . . . N

wherein p and q are two random solutions, and ϵ conforms to an uniform distribution in the interval [0,1].

After completing the global update and partial update, return to step 2 and the iterative operation is performed in sequence until the number of updates n reaches the preset maximum number of times N₀, that is, n≥N₀, and then go to step 4.

Step 4: The control variable x_(i) corresponding to Δs (x_(i))_(nmin) , the internal minimum deployment error value Δs(x_(i))_(nmin) and the corresponding trajectory are output. FIG. 6 shows the optimized trajectory of the target function for the minimum deployment error in this embodiment.

(6) A computer-readable/writable storage medium which stores a computer program.

When the computer program is executed, the computer program performs the following steps:

Step 1.: randomly generate values for N sets of initial control variables such as velocity, accelerations and orientations within a certain search range. Input all of N sets of xi to the sensor's underwater glide kinematics and dynamics models, and calculate the smallest distance between N actual deployment positions and target deployment positions. 3. The disclosure can improve the accuracy of prediction on the underwater three-dimensional trajectory of the marine sensor.

Step 2: input the selected N sets of initial control variables into the sensor's underwater glide kinematics model.

Step 3: calculate the smallest distance between N actual deployment positions and target deployment positions.

Step 4: determine whether the number of iterative operations is larger than the pre-set value, if yes, output the trajectory and the corresponding control variable as the optimal trajectory and optimal control variable. Otherwise, perform global random walk and local random walk on N sets of control variable, obtain N sets of control variable again, and return to step 2.

Step 5: save the optimal trajectory and optimal control variable such as initial velocity, acceleration and orientation output in step 4 into the computer-readable/writable storage medium.

(7) A trajectory optimization device for accurately deploying a marine sensor under water is provided, which includes the computer-readable/writable storage medium and a processor as described above, and the computer-readable/writable storage medium is able to store the computer program and the data output by the execution of the computer program, the processor is used for calling and processing computer programs stored in the computer-readable storage medium. Once the computer program is called, it starts to solve the sensor's underwater glide kinematics model, output the underwater trajectory of the sensor and the deployment distance error. Finally, the optimal trajectory and optimal control variable such as initial velocity, acceleration and orientation can be obtained by sequential global and local random updates for control variable such as the initial velocity, acceleration and orientation.

In summary, the trajectory optimization method for accurate underwater deployment of marine environment detection sensors provided by the disclosure consists of a set of underwater three-dimensional trajectory prediction models and trajectory optimization algorithms. The underwater three-dimensional trajectory prediction model includes the underwater glide kinematics model and dynamics model of the marine environment detection sensor. The underwater glide kinematics model mainly includes the underwater gliding motion equation of the marine environment detection sensor and the conversion relationship between the local coordinate system and the linear velocity and angular velocity of the global coordinate system obtained through Euler angles. The dynamics model refers to the analysis and description of the force on the marine environment detection sensor when performing gliding motion under water. The trajectory optimization algorithm refers to the use of modern optimization algorithms to find the optimal input conditions and optimal trajectories that meet the underwater deployment operations. The input of the trajectory optimization algorithm is the main initial control variables of the three-dimensional trajectory prediction model, including the initial linear velocity, angular velocity and Euler angle at the moment when the marine environment detection sensor is released. The output of the algorithm is the trajectory of underwater gliding motion of the marine environment detection sensor, the deployment position error and so on.

The underwater three-dimensional trajectory prediction model of the marine environment detection sensor and a trajectory optimization device provided by the disclosure can accurately describe and predict the underwater gliding motion of the marine environment detection sensor with six degrees of freedom. By taking into consideration the lateral force and displacement generated by the axial rotation or the asymmetric vortex, the accuracy of prediction on the underwater gliding trajectory of the marine environment detection sensor can be improved significantly. The trajectory optimization algorithm provided by the disclosure has flexible iteration steps, which can be large or small, free from the limitation of a single step length, so that the optimal solution search is efficient, and the trajectory can be optimized through multiple input variables.

In addition, the method for optimizing the underwater trajectory of a marine environment detection sensor and a trajectory optimization device provided by the disclosure can be applied to trajectory optimization under the condition of a wider range of deployment operations, including high speed, short range, etc., and can be used to provide a more accurate guidance for optimal design and implementation of underwater deploying operation of marine environment detection sensors. In the meantime, the underwater trajectory optimization method can be used for performing accurate deployment and arrangement for other underwater long-term observation type unpowered devices such as underwater seismic monitoring, acoustic detection sensors, and so on.

Those skilled in the art can easily understand that the above are only preferred embodiments of the disclosure and are not intended to limit the disclosure. Any modification, equivalent replacement and improvement, etc., made within the spirit and principle of the disclosure should be included in the protection scope of the disclosure. 

What is claimed is:
 1. A trajectory optimization method for accurately deploying marine sensors under water, comprising the following steps: step 1: according to a preset search range x_(min)≤x≤x_(max), N sets of values are randomly selected for initial control variables x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(z0), ϕ₀, ψ₀), denoted as x_(i), wherein, i=1,2 . . . , N, representing an i-th set, x_(min), x_(max) are a lower limit and an upper limit of the search range, U₀, V₀ and W₀ are components of an initial linear velocity in x-, y-, z-directions in an local coordinate system, ω_(x0), ω_(y0) and ω_(z0) are components of an initial angular velocity ω₀ in the x-, y-, z-directions in the local coordinate system, and (ϕ₀, θ₀, ψ₀) are initial Euler angles of the sensor; step 2: in an n-th iterative operation, all N sets of x_(i) are input into a sensor's underwater glide kinematics and dynamics model, and N actual deployment positions are calculated, then, a minimum deployment error Δs(x_(i))_(min) between a target deployment position and the actual deployment position of each set of x_(i) in the current iterative operation is calculated, thereafter, minimum distances corresponding to the N sets of x_(i) are compared with each other, a minimum value is taken as a global minimum deployment error Δs(x)_(nminmin) of this iterative operation, wherein: when n=1, Δs(x_(i))_(1min)=Δs(x_(i))₁, Δs(x)_(1minmin) is denoted as follows, and turn to step 3: Δs(x)_(1minmin)=Min{Δs(x _(i))₁ , Δs(x ₂)₁, . . . Δs(x _(N))₁}; when n>1, Δs(x_(i))_(nmin)=Min{Δs(x_(i))_(n), Δs(x_(i))_((n−1))min}, and then turn to step 3, under the circumstances, Δs(x)_(nminmin)=Min{Δs(x ₁)_(nmin), Δs(x ₂)_(nmin), . . . Δs(x _(N))_(nmin)} wherein a subscript n represents current iteration operation steps, Δs(x_(i))_(n), represents an error distance between the actual deployment position and the target deployment position in the n-th step of iterative operation in the i-th set, a subscript min is a mark for a minimum value in the set, and a subscript minmin is a mark for minimum value between the sets, Δs(x_(i))_(nmin) represents the internal minimum deployment error in the n-th step of iterative operation in the i-th set, and Δs (x)_(nminmin) represents a global minimum deployment error of a whole domain in the n-th step of iterative operation; step 3: determine whether n is less than a preset maximum number of iteration steps N₀, if yes, a global random walk is performed, N sets of x_(i), are updated, and then a local random walk is performed on the above updated N sets of x_(i) according to a preset probability, the N sets of x are obtained again, and return to step 2; otherwise, go to step 4; step 4: the control variable x_(i) corresponding to Δs(x)_(nminmin) and a corresponding trajectory are output as an optimal control variable and an optimal trajectory.
 2. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 1, wherein, in step 2, a kinematics model of the sensor's underwater gliding motion is as follows: $\overset{.}{U} = {{\frac{1}{m}\left\lbrack {{\left( {G - B} \right)\sin\theta} + F_{dx} + F_{ex}} \right\rbrack} - {W\omega_{y}} + {V\omega_{Z}} + {L_{G}\left( {\omega_{y}^{2} + \omega_{z}^{2}} \right)}}$ $\overset{.}{V} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\theta\sin\phi} + F_{dy} + F_{ly} + F_{ey}} \right\rbrack} - {U\omega_{z}} - {L_{G}{\overset{.}{\omega}}_{z}}}$ $\overset{.}{W} = {{\frac{1}{m}\left\lbrack {{{- \left( {G - B} \right)}\cos\theta\cos\phi} + F_{dz} + F_{lz} + F_{ez}} \right\rbrack} + {U\omega_{y}} + {L_{G}{\overset{.}{\omega}}_{y}}}$ ${\overset{.}{\omega}}_{y} = {\frac{1}{J_{y}}\left( {M_{ly} + M_{dy} + M_{ey} + M_{GBy} + {mL_{G}\overset{\overset{.}{\;}}{W}} - {mL_{G}U\omega_{y}}} \right)}$ ${\overset{.}{\omega}}_{z} = {\frac{1}{J_{z}}\left( {M_{lz} + M_{dz} + M_{ez} - {mL_{G}\overset{.}{V}} - {mL_{G}U\omega_{z}}} \right)}$ in the formula, {dot over (U)}, {dot over (V)}, and {dot over (W)} are acceleration components of the sensor in the x-, y-, z-directions in the local coordinate system; {dot over (ω)} _(y) and {dot over (ω)} _(z) are angular acceleration components of the sensor in the y-, z-directions in the local coordinate system; m is a mass of the sensor, U, V, and W are velocity components of the sensor in the x-, y-, z-directions in the local coordinate system, ω_(x), ω_(y) and ω_(z) are angular velocities of the sensor in the x-, y-, z-directions in the local coordinate system, ϕ, θ, and ψ are roll, pitch and yaw angles of the sensor in the local coordinate system, J_(y) and J_(z) are moment of inertia of the sensor with respect to the y-axis and z-axis, respectively, B is a buoyancy force of the sensor, M_(GBy) is a hydrostatic moment caused by a relative offset L_(G) between the buoyancy center B and a gravity center G, F_(e) and M_(e) are a hydrodynamic force and a moment of the sensor, F_(l) and M_(l) are a lifting force and a torque received by the sensor, F_(d) and M_(d) are a drag force and a torque received by the sensor, the subscripts x, y, z indicate that the corresponding parameter is a component on the corresponding coordinate axis; a relationship between the angular velocity and the Euler angles is obtained through a conversion rule for a local coordinate system and a global coordinate system: $\begin{bmatrix} \omega_{x} \\ \omega_{y} \\ \omega_{z} \end{bmatrix} = {\begin{bmatrix} {{- \overset{.}{\psi}}{\sin(\theta)}} \\ {\overset{.}{\psi}{\sin(\phi)}{\cos(\theta)}} \\ {\overset{.}{\psi}{\cos(\phi)}{\cos(\theta)}} \end{bmatrix} + \begin{bmatrix} 0 \\ {\overset{\overset{.}{\;}}{\theta}{\cos(\phi)}} \\ {{- \overset{\overset{.}{\;}}{\theta}}{\sin(\phi)}} \end{bmatrix} + \begin{bmatrix} \overset{.}{\phi} \\ 0 \\ 0 \end{bmatrix}}$ wherein {dot over (ψ)}, {dot over (θ)}, {dot over (ϕ)} are change rates of a yaw angle ψ, a pitch angle θ and a roll angle ϕ that varies with time.
 3. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein the hydrodynamic forces in the y and z directions are calculated by the following expressions: F _(ey) =−{dot over (V)}m ₂₂ −{dot over (ω)} _(y) m ₂₄ −UVM ₂₂(x _(t))+Uω _(z) x _(t) M ₂₂ (x _(t)) F _(ez) =−{dot over (W)}m ₃₃ −{dot over (ω)} _(z) m ₃₅ −UVM ₃₃(x _(t))+Uω _(y) x _(t) M ₃₃ (x _(t)) hydrodynamic torques in the y and z directions are calculated according to the following formula: M _(ey) =−{dot over (W)}m ₃₅−{dot over (ω)} _(y) m ₅₅ −UW[m _(mm) +x _(t) mM ₃₃(x _(t))]+Uω _(y) [m ₃₅−(x _(t))² M ₃₃(x _(t))] M _(ez) =−{dot over (V)}m ₂₆ −{dot over (ω)} _(x) m ₆₆ −UV[m ₂₂ +x _(t) M ₂₂(x _(t))]+Uω _(z) [m ₂₆−(x _(t))² M ₂₂(x _(t))] in the formula, m₂₂, m₂₄, m₂₆, m₃₃, m₃₅, m₅₅ and m₆₆ are three-dimensional additional mass coefficients of the sensor, and M₂₂ (x_(t)) and M₃₃ (x_(t)) respectively represent two-dimensional additional mass coefficients of an effective trailing edge in the y and z directions, x_(t) is a longitudinal position of an effective trailing edge; for a sensor with a tip, M₃₃ (x_(t))=0 M₂₂ (x_(t))=0, for a marine sensor with two blunt ends, M₃₃ (X_(t))≠0, M₂₂ (X_(t))≠0.
 4. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein the drag force and the torque are generated by a relative translational movement of the sensor in contact with water, and an axial drag force F_(dx) is calculated as follows: $F_{dx} = \left\{ \begin{matrix} {{{{- C_{1}}\pi D\sqrt{v\rho^{2}L}U\sqrt{|U|}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{{Laminar}\mspace{14mu}{flow}^{''}}} \\ {{{{- C_{2}}\frac{1}{2}\rho\pi DLU^{2}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{{Transition}\mspace{14mu}{flow}^{''}}} \\ {{{C_{3}\frac{1}{2}\rho\pi DLU^{2}} - {\frac{1}{8}\rho\pi C_{dt}D^{2}U{U}}},} & {\;^{''}{Turbulence}^{''}} \end{matrix} \right.$ wherein D is a cross-sectional diameter of the sensor, v is a kinematic viscosity of the water, ρ is a density of the water, C_(dt) is a tangential shape drag coefficient of the sensor, C₁, C₂, and C₃ are shape drag coefficients of the sensor when there are laminar flow, transition flow and turbulence, L is a length of the sensor; lateral drag forces F_(dy) and F_(dz) in the y and z directions are as calculated follows: F _(dy)=½ρC _(dy) D ∫ _(−0.5L) ^(0.5L)−(V+ω _(z) x|dx F _(dz)=½ρC _(dz) D ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)|V+ω _(z) x|dx resistance moments M_(dy) and M_(dz) caused by the lateral drag force are calculated as follows: M _(dy)=½ρC _(dz) D ∫ _(−0.5L) ^(0.5L) x(w−ω _(y) x)|W−ω _(y) x|dx M _(dz)=−½ρC _(dy) D ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)|V+ω _(z) x|dx wherein C_(dy) and C_(dz) are respectively the drag coefficients in the y and z directions.
 5. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein if the marine sensor itself does not rotate around the axis, the lifting force F_(l) and the moment M_(l) are the lateral forces that are mainly caused by the asymmetry of a wake vortex distribution during the flow of water laterally crossing the sensor; if the sensor itself rotates, the water flow circulates laterally and generates pressure difference between two sides of the sensor, which forms the lifting force; the lifting force that is generated when the sensor does not rotate in the axial direction is caused by the asymmetry of the wake vortex distribution, and the lifting forces are calculated as follows: F _(ly)=½ρD ∫ _(−0.5L) ^(0.5L) (W−ω _(y) x)C _(l)(x)dx F _(lz)=−½ρD ∫ _(−0.5L) ^(0.5L) (V+ω _(z) x)C _(l)(x)dx M _(ly)=½ρD ∫ _(−0.5L) ^(0.5L) x(V+ω _(z) x)C _(l)(x)dx M _(lz)=½ρD ∫ _(−0.5L) ^(0.5L) x(W−ω _(y) x)C _(l)(x)dx wherein the lifting force coefficient C_(l)(x) is calculated as follows: ${C_{l}(x)} = \left\{ {\begin{matrix} {{KC_{La}{\sin\left\lbrack {\pi\;{{St}\left( {t^{\prime} - 8} \right)}} \right\rbrack}},} & {t^{\prime} > 8} \\ {0,} & {t^{\prime} \leq 8} \end{matrix},} \right.$ wherein St=f_(v)D/V_(a) is the Strouhal number, f_(v) is a vortex shedding frequency, V_(a) is a cross-flow velocity, K=±1 reflects a random behavior of the lifting direction during gliding, C_(La) is a lift constant, and the time coefficient t′ is calculated as follows: $\begin{matrix} {t^{\prime} = {\frac{s(x)}{R} - {{0.3}51}}} & \; \end{matrix}$ wherein s(x) is a displacement in the y direction when the water flows around a cross section of the sensor, and R is a radius of the cross section of the sensor; the lifting force and moment during axial rotation are calculated as follows: F _(ly)½ρDW ∫ _(−0.5L) ^(0.5L) C′ _(l) (x)(W−ω _(y) x)dx F _(lz)=−½ρDV ∫ _(−0.5L) ^(0.5L) C′ _(l) (x)(V+ω _(z) x)dx M _(ly)=½ρDV ∫ _(−0.5L) ^(0.5L) C′ _(l) (x)(V+ω _(z) x)xdx M _(lz)=½ρDC′ _(l) (x)W ∫_(−0.5L) ^(0.5L) (W−ω _(y) x)xdx wherein the lifting force coefficient C′_(l) (x) is calculated as follows: $\begin{matrix} {{C_{l}^{\prime}(x)} = \left\{ {\begin{matrix} {\frac{\omega_{x}D}{V_{a}},} & {\frac{\omega_{x}D}{V_{a}} \leq 8} \\ {{8 + {0.12\left( {\frac{\omega_{x}D}{v_{a}} - 8} \right)}},} & {\frac{\omega_{x}D}{V_{a}} > 8} \end{matrix}.} \right.} & \; \end{matrix}$
 6. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 1, wherein the sensor in the local body coordinate system is an axisymmetric slender body, and the change of an initial roll angle does not affect the movement of the slender body, therefore, the initial roll angle and an initial immersion depth are as follows: ϕ₀=0 d=−L/2sin(θ₀) thus, the initial control variable is simplified as x=(U₀, V₀, W₀, ω_(x0), ω_(z0), θ₀, ψ₀), and L represents the length of the sensor.
 7. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 1, wherein in step 3, a global random walk is performed through Lévy flight, and global update is performed for N sets of control parameters x_(i), as shown in the following formula: x _(i) ^(n+1) =x _(i) ^(n) +a ⊗ Lévy(s, λ) in the formula, a=a₀ (x_(j) ^(t)−x_(i) ^(t)), a₀ is a preset constant, x_(i) represents the j-th set of control parameters x that is randomly generated, s represents the step length of the iterative operation, ⊗ stands for item-based multiplication, and λ is an exponent of a Lévy distribution.
 8. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 1, wherein after the global update in step 3, each set of control parameters x_(i) has a probability of p_(a) ∈(0,1) to be locally update, and the update method of the local random walk is shown in the following formula: x _(i) ^(n+1) =x _(i) ^(n+1)+ϵ·(x _(p) ^(n+1) −x _(q) ^(n+1)), p, q=1, 2, . . . N wherein x_(p) and x_(q) are two random solutions, ϵ conforms to an uniform distribution of an interval [0,1], and the superscript n+1 is the number of iterations in the next step.
 9. A computer-readable/writable storage medium storing a computer program, wherein when the computer program is executed, the computer program performs the following steps: step 1: randomly generating values for N sets of initial control variables such as velocity, accelerations and orientations within a certain search range; inputting all of N sets of xi to the sensor's underwater glide kinematics and dynamics models, and calculating a smallest distance between N actual deployment positions and target deployment positions; step 2: inputting the selected N sets of initial control variables into the sensor's underwater glide kinematics model; step 3: calculating the smallest distance between N actual deployment positions and target deployment positions; step 4: determining whether a number of iterative operations is larger than a pre-set value, if yes, outputting the trajectory and the corresponding control variable as an optimal trajectory and an optimal control variable, otherwise, performing global random walk and local random walk on the N sets of control variable, obtaining N sets of control variable again, and returning to step 2; step 5: save the optimal trajectory and the optimal control variable such as initial velocity, acceleration and orientation output in step 4 into the computer-readable/writable storage medium.
 10. A trajectory optimization device for accurately deploying a marine sensor under water, comprising: a computer-readable/writable storage medium of claim 9 and a processor, wherein the computer-readable/writable storage medium is able to store the computer program and the data output by an execution of the computer program, the processor is used for calling and processing computer programs stored in the computer-readable storage medium, once the computer program is called, it starts to solve the sensor's underwater glide kinematics model, output the underwater trajectory of the sensor and the deployment distance error, finally, the optimal trajectory and optimal control variable such as initial velocity, acceleration and orientation can be obtained by sequential global and local random updates for control variable such as the initial velocity, acceleration and orientation.
 11. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein the sensor in the local body coordinate system is an axisymmetric slender body, and the change of an initial roll angle does not affect the movement of the slender body, therefore, the initial roll angle and an initial immersion depth are as follows: ϕ₀=0 d=−L/2sin(θ₀) thus, the initial control variable is simplified as x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(x0), θ₀, ψ₀), and L represents the length of the sensor.
 12. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein in step 3, a global random walk is performed through Lévy flight, and global update is performed for N sets of control parameters x_(i), as shown in the following formula: x _(i) ^(n+1) =x ₁ ^(n) +a ⊗ Lévy(s, λ) in the formula, a=a₀ (x_(j) ^(t)−x_(i) ^(t)), a₀ is a preset constant, x_(j) represents the j-th set of control parameters x that is randomly generated, s represents the step length of the iterative operation, ⊗ stands for item-based multiplication, and λ is an exponent of a Lévy distribution.
 13. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 2, wherein after the global update in step 3, each set of control parameters x_(i) has a probability of p_(a) ∈(0,1) to be locally update, and the update method of the local random walk is shown in the following formula: x _(i) ^(n+1) =x _(i) ^(n+1)+ϵ·(x _(p) ^(n+1) −x _(q) ^(n+1), p, q=1, 2, . . . N wherein x_(p) and x_(q) are two random solutions, ϵ conforms to an uniform distribution of an interval [0,1], and the superscript n+1 is the number of iterations in the next step.
 14. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, when the computer program is executed by a processor, the trajectory optimization method as claimed in claim 2 is implemented.
 15. A trajectory optimization device for accurately deploying a marine sensor under water, comprising the computer-readable storage medium and the processor as claimed in claim 14, and the processor is configured for calling and processing computer programs stored in the computer-readable storage medium.
 16. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 3, wherein the sensor in the local body coordinate system is an axisymmetric slender body, and the change of an initial roll angle does not affect the movement of the slender body, therefore, the initial roll angle and an initial immersion depth are as follows: ϕ₀=0 d=−L/2sin(θ₀) thus, the initial control variable is simplified as x=(U₀, V₀, W₀, ω_(x0), ω_(y0), ω_(z0), θ₀, ψ₀), and L represents the length of the sensor.
 17. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 3, wherein in step 3, a global random walk is performed through Lévy flight, and global update is performed for N sets of control parameters x_(i), as shown in the following formula: x ₁ ^(n+1) =x _(i) ^(n) +a ⊗ Lévy (s, λ) in the formula, a=a₀(x _(j) ^(t) −x _(i) ^(t)), a₀, is a preset constant, x_(j) represents the j-th set of control parameters x that is randomly generated, s represents the step length of the iterative operation, ⊗ stands for item-based multiplication, and λ is an exponent of a Lévy distribution.
 18. The trajectory optimization method for accurately deploying marine sensors underwater according to claim 3, wherein after the global update in step 3, each set of control parameters x_(i) has a probability of p_(a) ∈(0,1) to be locally update, and the update method of the local random walk is shown in the following formula: x _(i) ^(n+1) =x _(i) ^(n+1)+ϵ·(x _(p) ^(n+1) −x _(q) ^(n+1)), p, q=1, 2, . . . , N wherein x_(p) and x_(q) are two random solutions, c conforms to an uniform distribution of an interval [0,1], and the superscript n+1 is the number of iterations in the next step.
 19. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, when the computer program is executed by a processor, the trajectory optimization method as claimed in claim 3 is implemented.
 20. A trajectory optimization device for accurately deploying a marine sensor under water, comprising the computer-readable storage medium and the processor as claimed in claim 19, and the processor is configured for calling and processing computer programs stored in the computer-readable storage medium. 